/* jerasure-2.cpp
 * James S. Plank

 Jerasure - A C/C++ Library for a Variety of Reed-Solomon and RAID-6 Erasure Coding Techniques

 Revision 2.0
 May 24, 2011

 James S. Plank
 Department of Electrical Engineering and Computer Science
 University of Tennessee
 Knoxville, TN 37996
 plank@cs.utk.edu

 Copyright (c) 2011, James S. Plank
 All rights reserved.

 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions
 are met:

 - Redistributions of source code must retain the above copyright
 notice, this list of conditions and the following disclaimer.

 - Redistributions in binary form must reproduce the above copyright
 notice, this list of conditions and the following disclaimer in
 the documentation and/or other materials provided with the
 distribution.

 - Neither the name of the University of Tennessee nor the names of its
 contributors may be used to endorse or promote products derived
 from this software without specific prior written permission.

 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
 OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
 AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY
 WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 POSSIBILITY OF SUCH DAMAGE.

 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "galois-2.h"
#include "jerasure-2.h"

#define talloc(type, num) (type *) malloc(sizeof(type)*(num))

static double jerasure_total_xor_bytes = 0;
static double jerasure_total_gf_bytes = 0;
static double jerasure_total_memcpy_bytes = 0;

/**** JER_Bitmatrix methods: ****/

//constructors:
JER_Bitmatrix::JER_Bitmatrix(int superr, int superc){
  int i;

  W = 8;
  SUPERC = superc;
  /* Matrix in row major form. */
  Elts.resize(superr*W);

  for( i = 0; i < Elts.size(); i++ )
  {
    /* Each row is a vector of characters. Divide by sizeof( char ),
       which should always be 1 but not sure about portability, times
       8, number of bits in byte. */
    Elts[i].resize( ((superc*W) - 1) / (8*sizeof(char)) + 1, 0 );
  }
}

JER_Bitmatrix::JER_Bitmatrix(int superr, int superc, int w){
  int i;

  W = w;
  SUPERC = superc;
  /* Matrix in row major form. */
  Elts.resize(superr*W);

  for( i = 0; i < Elts.size(); i++ )
  {
    /* Each row is a vector of characters. Divide by sizeof( char ),
       which should always be 1 but not sure about portability, times
       8, number of bits in byte. */
    Elts[i].resize( ((superc*W) - 1) / (8*sizeof(char)) + 1, 0 );
  }
}

JER_Bitmatrix::JER_Bitmatrix(JER_Matrix &jm){
  int superr, superc;
  int i, j, k, x, l;
  unsigned int elt;

  W = jm.W;

  /* If superr <= 0, matrix will not be of valid size. But if it is
     okay to use an invalid JER_Matrix, it should be okay to use an
     invalid JER_Bitmatrix. Checks prevent code from crashing. */
  superr = jm.Elts.size();
  if( superr > 0 )
  {	
    superc = jm.Elts[0].size();
    SUPERC = superc;
  }
  else
    return;

  /* Set up vectors to appropriate size. */
  try{ 
	  Elts.resize( superr * W ); 
  for( i = 0; i < Elts.size(); i++ )
    Elts[i].resize( 
        ((superc*W) - 1) / (8*sizeof(char)) + 1, 0 );
  }
  catch(...){ return; }

  /* Copied more or less from original. */
  for( i = 0; i < superr; i++ )
  {
    for( j = 0; j < superc; j++ )
    {
      /* Take the value at matrix[i][j] and store it as it bitmatrix
         specification. */
      elt = jm.Elts[i][j];
      for( x = 0; x < W; x++ )
      {
        /* Set the bit in the xth row of the ith super row. */
        for( l = 0; l < W; l++ )
        {
          /* Set the lth bits of the jth super col to the bit rep
             of elt. */
          Set( i*W+l, j*W+x, ((elt & (1 << l)) ? 1 : 0) );
        }
        elt = (unsigned int) galois_single_multiply( elt, 2, W );
      }
    }
  }
}

JER_Bitmatrix::JER_Bitmatrix( JER_Bitmatrix &jbm )
{
  int i;
  W = jbm.W;
  SUPERC = jbm.SUPERC;

  Elts.resize( jbm.Elts.size() );

  for( i = 0; i < Elts.size(); i++ )
  {
    Elts[i] = jbm.Elts[i];
  }
}

/* Sets the bit at (r,c) to val. */
void JER_Bitmatrix::Set(int r, int c, int val){
  int index, offset;
  char mask;

  mask = -1;

  index = c / ( 8 * sizeof(char) );
  offset = c % ( 8 * sizeof(char) );

  if( val )
    Elts[r][index] = Elts[r][index] | (1 << offset);
  else
    Elts[r][index] = Elts[r][index] & ( mask ^ ( 1 << offset ) );
}

/* Changed from void. Should return 0 or 1. */
int JER_Bitmatrix::Isset(int r, int c){
  int index, offset;

  index = c / ( 8 * sizeof(char) );
  offset = c % ( 8 * sizeof(char) );

  return ( Elts[r][index] & ( 1 << offset ) ) >> offset;
}

void JER_Bitmatrix::CreateIdentity()
{
  int i, j;

  if( Elts.size() == SUPERC * W )
  {
    for( i = 0 ; i < Elts.size(); i++ )
    {
      Elts[i].clear();
      Elts[i].resize( ( (SUPERC * W - 1) / (8 * sizeof(char)) ) + 1 );
      Set( i, i, 1);
    }
  }
}

int JER_Bitmatrix::Invertible(){
  JER_Bitmatrix* jbm;
  int rows, cols;
  int i, j, k;
  
  jbm = new JER_Bitmatrix( *((JER_Bitmatrix*)this) );
  if( jbm == NULL )
  {
    return -1;
  }

  rows = Elts.size();
  cols = SUPERC * W;

  if( rows != cols )
	  return 0;

  /* First -- convert int upper triangular */
  for( i = 0; i < cols; i++ )
  {
    if( !jbm->Isset( i, i ) )
    {
      for( j = i + 1; j < rows && !jbm->Isset(j,i); j++ );
      if( j == rows ) 
      {
        delete jbm;
        return 0;
      }
      jbm->Elts[i].swap(jbm->Elts[j]);
    }

    for( j = i + 1; j < rows; j++ )
    {
      if( jbm->Isset( j, i ) )
      {
        for( k = 0; k < Elts[0].size(); k++ )
        {
          jbm->Elts[j][k] ^= jbm->Elts[i][k];
        }
      }
    }
  }

  delete jbm;
  return 1;
}

int JER_Bitmatrix::Inverse(){
  JER_Bitmatrix *jbm;
  int rows, cols;
  int i, j, k;

  jbm = new JER_Bitmatrix( *((JER_Bitmatrix*)this) );
  if( jbm == NULL )
  {
    return -1;
  }
  CreateIdentity();

  rows = Elts.size();
  cols = SUPERC * W;

  if( rows != cols )
	  return -1;

  /* First -- convert int upper triangular */
  for( i = 0; i < cols; i++ )
  {
    if( !jbm->Isset( i, i ) )
    {
      for( j = i + 1; j < rows && !jbm->Isset(j,i); j++ );
      if( j == rows ) 
      {
        delete jbm;
        return -1;
      }
      jbm->Elts[i].swap(jbm->Elts[j]);
      Elts[i].swap(Elts[j]);
    }

    for( j = i + 1; j < rows; j++ )
    {
      if( jbm->Isset( j, i ) )
      {
        for( k = 0; k < Elts[0].size(); k++ )
        {
          jbm->Elts[j][k] ^= jbm->Elts[i][k];
          Elts[j][k] ^= Elts[i][k];
        }
      }
    }
  }

  /* Now the matrix is upper triangular. Start at the top and multiply down. */

  for( i = rows - 1; i >= 0; i-- )
  {
    for( j = 0; j < i; j++ )
    {
      if( jbm->Isset(j,i) )
      {
        for( k = 0; k < Elts[0].size(); k++ )
        {
          jbm->Elts[j][k] ^= jbm->Elts[i][k];
          Elts[j][k] ^= Elts[i][k];
        }
      }
    }
  }
  delete jbm;
  return 0;
}

void JER_Bitmatrix::Print(){
  int i, j;

  //loop through all rows
  for (i = 0; i < Elts.size(); i++) {

    //insert a new line separating super-rows
    if (i != 0 && i%W == 0){
      printf("\n");
    }

    //loop through every column
    for (j = 0; j < SUPERC*W; j++) {

      //print spaces between super columns
      if (j != 0 && j%W == 0) printf(" ");

      printf("%d", Isset(i,j));

    }
    printf("\n");
  }

}
string JER_Bitmatrix::String(){
  string s;
  int r,c;

  for (r = 0; r < Elts.size(); r++) {
    for(c=0;c<SUPERC*W;c++){
      if(Isset(r,c)){
        s+="1";
      }else{
        s+="0";
      }
    }
  }

  return s;

}

/**** JER_Matrix methods: ****/

JER_Matrix::JER_Matrix( int r, int c, int wordsize )
{
  int i;

  W = wordsize;
  Elts.resize( r );
  for( i = 0; i < r; i++ )
    Elts[i].resize( c, 0 );
}

JER_Matrix::JER_Matrix( JER_Matrix &jm )
{
  int i;

  W = jm.W;
  Elts.resize( jm.Elts.size() );

  for( i = 0; i < Elts.size(); i++ )
    Elts[i] = jm.Elts[i];
}

void JER_Matrix::CreateIdentity()
{
  int i;
  int rows;

  rows = Elts.size();
  if( rows > 0 && rows == Elts[0].size() )
  {
    for( i = 0; i < rows; i++ )
    {
      Elts[i].clear();
      Elts[i].resize( rows, 0 );
      Elts[i][i] = 1;
    }
  }
}

int JER_Matrix::Invertible()
{
  JER_Matrix* jm;
  int cols, rows;
  int i, j, k, tmp;

  rows = cols = Elts.size();
  if( rows <= 0 || Elts[0].size() != rows )
    return 0;

  jm = new JER_Matrix( *((JER_Matrix*)this) );
  if( jm == NULL )
    return 0;

  for( i = 0; i < rows; i++ )
  {
    if( jm->Elts[i][i] == 0 )
    {
      for( j = i + 1; j < rows && jm->Elts[j][i] == 0; j++ );
      if( j == rows ){ delete jm;	return 0; }
      jm->Elts[i].swap( jm->Elts[j] );
    }
    tmp = jm->Elts[i][i];
    if( tmp != 1 )
    {
      tmp = galois_single_divide( 1, tmp, W );
      for( j = 0; j < cols; j++ )
      {
        jm->Elts[i][j] = galois_single_multiply( jm->Elts[i][j], tmp, W );
      }
    }
    for( j = i + 1; j < rows; j++ )
    {
		/* Every row below i is multiplied by row i times [j,i] so that
		   all entries in column i below row i = 0. */
      if( jm->Elts[j][i] != 0 )
      {
        if( jm->Elts[j][i] == 1 )
        {
          for( k = 0; k < cols; k++ )
          {
            jm->Elts[j][k] ^= jm->Elts[i][k];
          }
        }
        else
        {
          tmp = jm->Elts[j][i];
          for( k = 0; k < cols; k++ )
          {
            jm->Elts[j][k] ^= galois_single_multiply( 
                tmp, jm->Elts[i][k], W );
          }
        }
      }
    }
  }
  delete jm;
  return 1;
}

int JER_Matrix::Inverse()
{
  JER_Matrix* jm;
  int cols, rows;
  int i, j, k, tmp;

  rows = cols = Elts.size();
  if( rows <= 0 || Elts[0].size() != rows )
    return -1;

  jm = new JER_Matrix( *((JER_Matrix*)this) );
  if( jm == NULL )
    return -1;
  CreateIdentity();

  for( i = 0; i < rows; i++ )
  {
	  /* If there is a 0 on the diagonal on the ith row */
    if( jm->Elts[i][i] == 0 )
    {
		/* Go through each subsequent row looking for one with a non-zero
		   value in the ith column. When this row is moved to the ith row,
		   there will no longer be a 0 on the diagonal on the ith row. */
      for( j = i + 1; j < rows && jm->Elts[j][i] == 0; j++ );
	  /* If j == rows, there was no row with a non-zero value in the ith
		 column below row i, therefore matrix is not invertible. */
      if( j == rows ){ delete jm;	return -1; }
	  /* Otherwise swap the ith and jth rows. */
      jm->Elts[i].swap( jm->Elts[j] );
      Elts[i].swap( Elts[j] );
    }

	/* If the value on the diagonal is not 1, it must be reduced to 1 by 
	   dividing the row by the value on the diagonal. */
    tmp = jm->Elts[i][i];
    if( tmp != 1 )
    {
      tmp = galois_single_divide( 1, tmp, W );
      for( j = 0; j < cols; j++ )
      {
        jm->Elts[i][j] = galois_single_multiply( jm->Elts[i][j], tmp, W );
        Elts[i][j] = galois_single_multiply( Elts[i][j], tmp, W );
      }
    }

	/* jm->Elts[i][i] == 1 now. */
    for( j = i + 1; j < rows; j++ )
    {
		/* Every row below i is multiplied by row i times [j,i] so that
		   all entries in column i below row i = 0. */
      if( jm->Elts[j][i] != 0 )
      {
        if( jm->Elts[j][i] == 1 )
        {
          for( k = 0; k < cols; k++ )
          {
            jm->Elts[j][k] ^= jm->Elts[i][k];
            Elts[j][k] ^= Elts[i][k];
          }
        }
        else
        {
          tmp = jm->Elts[j][i];
          for( k = 0; k < cols; k++ )
          {
            jm->Elts[j][k] ^= galois_single_multiply( 
                tmp, jm->Elts[i][k], W );
            Elts[j][k] ^= galois_single_multiply( 
                tmp, Elts[i][k], W );
          }
        }
      }
    }
  }

  /* Now the matrix is upper triangular.  Start at the top and multiply down  */

  for( i = rows - 1; i >= 0; i-- )
  {
    for( j = 0; j < i; j++ )
    {
      if( jm->Elts[j][i] != 0 )
      {
        tmp = jm->Elts[j][i];
        for( k = 0; k < cols; k++ )
        {
          Elts[j][k] ^= galois_single_multiply( tmp, Elts[i][k], W );
        }
      }
    }
  }
  delete jm;
  return 0;
}

void JER_Matrix::Print(){
  int i, j;
  int fw;
  char s[30];
  unsigned int w2;
  int rows = Elts.size();
  int cols;

  if(rows){
    cols = Elts[0].size();
  }

  if (W == 32) {
    fw = 10;
  } else {
    w2 = (1 << W);
    sprintf(s, "%u", w2-1);
    fw = strlen(s);
  }

  for (i = 0; i < rows; i++) {
    for (j = 0; j < cols; j++) {
      if (j != 0) printf(" ");

      printf("%*u", fw, Elts[i][j]);
    }
    printf("\n");
  }

}

string JER_Matrix::String(){
  int r,c;
  string s="";
  int rows;
  int cols;
  char sn[11];

  rows = Elts.size();
  if(rows)
    cols = Elts[0].size();

  for(r=0;r<rows;r++){
    for(c=0;c<cols;c++){
      sprintf(sn,"%d",Elts[r][c]);
      s+= sn;
      s+=" ";
    }
  }

  //remove the last space
  if( s != "" )
    s.erase(s.end()-1);

  return s;
}

/***** Start schedule functions *****/

JER_Schedule *JER_BM_To_Schedule_Dumb(JER_Bitmatrix *jbm){

  JER_Schedule * sched;
  int r,c;
  int rows,cols;
  int optodo,op;
  IVec v5; //5 element vector

  rows = jbm->Elts.size();
  if(!rows) return NULL;

  cols = jbm->W*jbm->SUPERC;
  if(!cols) return NULL;

  sched = new JER_Schedule();
  if( sched == NULL )
	  return NULL;

  op = 0;
  for(r=0;r<rows;r++){
    optodo = 0;
    for(c=0;c<cols;c++){
      /* If jbm[r][c] == 1, xor (copy if its the first time) device c/w 
         packet c%w into device k + (r/w) packet r%w. */
      if(jbm->Isset(r,c)){
        sched->Elts.push_back(v5);
        sched->Elts[op].resize(5);
        sched->Elts[op][4] = optodo;
        sched->Elts[op][0] = c/jbm->W;
        sched->Elts[op][1] = c%jbm->W;
        sched->Elts[op][2] = jbm->SUPERC+r/jbm->W;
        sched->Elts[op][3] = r%jbm->W;
        optodo = 1;
        op++;
      }
    }
  }

  /* last schedules' 0th element (source disk) no longer needs to be -1, 
     because with vectors we can just use size() */
  return sched;
}

JER_Schedule *JER_BM_To_Schedule_CSHR(JER_Bitmatrix *jbm){

  JER_Schedule * sched;
  int r,c;
  int rows,cols;
  int optodo;
  IVec v5; //5 element vector
  int *diff, *from, *flink, *blink;
  int i, j;
  int no,row;
  int op;
  int bestrow, bestdiff, top;

  rows = jbm->Elts.size();
  if(!rows) return NULL;

  cols = jbm->W*jbm->SUPERC;
  if(!cols) return NULL;

  sched = new JER_Schedule();
  if( sched == NULL )
	  return NULL;

  diff = talloc(int, rows);
  from = talloc(int, rows);
  flink = talloc(int, rows);
  blink = talloc(int, rows);
  if( diff == NULL || from == NULL || flink == NULL || blink == NULL )
	  return NULL;

  bestdiff = cols+1;
  top = 0;
  op =0;
  for (i = 0; i < rows; i++) {
    no = 0;
    /* Count number of 1s in row i. Store in diff[i]. */
    for (j = 0; j < cols; j++) {
      no += jbm->Isset(i,j);
    }
    diff[i] = no;
    from[i] = -1;
    flink[i] = i+1;
    blink[i] = i-1;
    /* Find row with least number of 1s. */
    if (no < bestdiff) {
      bestdiff = no;
      bestrow = i;
    }
  }

  flink[rows-1] = -1;

  while (top != -1) {
    row = bestrow;

    /* Pull row out of double-linked list so that it may be processed. */
    if (blink[row] == -1) {
      top = flink[row];
      if (top != -1) blink[top] = -1;
    } else {
      flink[blink[row]] = flink[row];
      if (flink[row] != -1) {
        blink[flink[row]] = blink[row];
      }
    }


    /* Row does not benefit from having dependencies, simply generate it
       dumb. */
    if (from[row] == -1) {
      optodo = 0;
      for (j = 0; j < cols; j++) {
        if (jbm->Isset(row,j)) {
          sched->Elts.push_back(v5);
          sched->Elts[op].resize(5);
          sched->Elts[op][4] = optodo;
          sched->Elts[op][0] = j/jbm->W;
          sched->Elts[op][1] = j%jbm->W;
          sched->Elts[op][2] = jbm->SUPERC+row/jbm->W;
          sched->Elts[op][3] = row%jbm->W;
          op++;
          optodo = 1;
        }
      }
    }
    /* Otherwise generate it starting from its dependency. 
       The first operation is a copy of from[row]; each subsequent action
       is xor with individual bits. */
    else {
      sched->Elts.push_back(v5);
      sched->Elts[op].resize(5);
      sched->Elts[op][4] = 0;
      sched->Elts[op][0] = jbm->SUPERC + from[row]/jbm->W;
      sched->Elts[op][1] = from[row]%jbm->W;
      sched->Elts[op][2] = jbm->SUPERC+row/jbm->W;
      sched->Elts[op][3] = row%jbm->W;
      op++;

      for (j = 0; j < cols; j++) {
        if (jbm->Isset(row,j) ^ jbm->Isset(from[row],j)) {
          sched->Elts.push_back(v5);
          sched->Elts[op].resize(5);
          sched->Elts[op][4] = 1;
          sched->Elts[op][0] = j/jbm->W;
          sched->Elts[op][1] = j%jbm->W;
          sched->Elts[op][2] = jbm->SUPERC+row/jbm->W;
          sched->Elts[op][3] = row%jbm->W;
          op++;
          optodo = 1;
        }
      }
    }
    bestdiff = cols+1;

    /* Look for next smallest row. Before checking diff[i] for min, first
       see if diff[i] reduced by xoring the recently processed row with 
       row i. */
    for (i = top; i != -1; i = flink[i]) {
      /* Operations for using previous entry = 1 + number of 1s in xor of
         two rows. */
      no = 1;
      for (j = 0; j < cols; j++) no += (jbm->Isset(row,j) ^ jbm->Isset(i,j));
      if (no < diff[i]) {
        from[i] = row;
        diff[i] = no;
      }
      if (diff[i] < bestdiff) {
        bestdiff = diff[i];
        bestrow = i;
      }
    }
  }

  free(from);
  free(diff);
  free(blink);
  free(flink);

  return sched;
}

void JER_Do_Scheduled_Operations(vector <unsigned char *> &ptrs, JER_Schedule *sched,int packetsize){
  unsigned char *sptr;
  unsigned char *dptr;
  int op;
  int n0,n1,n2,n3;

  for (op = 0; op< sched->Elts.size(); op++) {
    n0 = (sched->Elts)[op][0];
    n1 = (sched->Elts)[op][1];
    n2 = (sched->Elts)[op][2];
    n3 = (sched->Elts)[op][3];

    sptr = ptrs[n0] + n1*packetsize;
    dptr = ptrs[n2] + n3*packetsize;
    if ((sched->Elts)[op][4]) {
      galois_region_xor(sptr, dptr, dptr, packetsize);
      jerasure_total_xor_bytes += packetsize;
    } else {
      memcpy(dptr, sptr, packetsize);
      jerasure_total_memcpy_bytes += packetsize;
    }
  }
}

void JER_Schedule_Encode(JER_Slices *slices, JER_Schedule *sched){
  vector <unsigned char*> ptr_copy;
  int i, j, tdone;
  int k,m;

  k = slices->K;
  m = slices->M;

  try{ ptr_copy.resize(k+m); }
  catch(...) { return; }
  for (i = 0; i < k+m; i++) ptr_copy[i] = slices->ptrs[i];
  for (tdone = 0; tdone < slices->size; tdone += slices->PS*slices->W) {
    JER_Do_Scheduled_Operations(ptr_copy, sched,slices->PS);
    for (i = 0; i < k+m; i++) ptr_copy[i] += (slices->PS*slices->W);
  }
}

static vector <unsigned char *> *JER_Set_Up_Ptrs_For_Scheduled_Decoding(JER_Slices * slices, vector <int> & erasures)
{
  int ddf, cdf;
  vector <int> erased;
  vector <unsigned char *> *ptrs;
  int i, j, x;
  int k = slices->K;
  int m = slices->M;

  ddf = 0;
  cdf = 0;
  //count how many data and coding disks have failed
  for (i = 0; i<erasures.size(); i++) {
    if (erasures[i] < k) ddf++; else cdf++;
  }

  JER_Erasures_To_Erased(k, m, erasures, erased);
  if (erased.size() != k+m) return NULL;

  /* Set up ptrs.  It will be as follows:

     - If data drive i has not failed, then ptrs[i] = data_ptrs[i].
     - If data drive i has failed, then ptrs[i] = coding_ptrs[j], where j is the 
     lowest unused non-failed coding drive.
     - Elements k to k+ddf-1 are data_ptrs[] of the failed data drives.
     - Elements k+ddf to k+ddf+cdf-1 are coding_ptrs[] of the failed data drives.

     The array row_ids contains the ids of ptrs.
     The array ind_to_row_ids contains the row_id of drive i.

     However, we're going to set row_ids and ind_to_row in a different procedure.
   */

  ptrs = new vector <unsigned char *>;
  try{ ptrs->resize(k+m);}
  catch(...){
    delete(ptrs);
    return NULL;
  }

  j = k;
  x = k;
  for (i = 0; i < k; i++) {
    if (erased[i] == 0) {
      (*ptrs)[i] = slices->ptrs[i];
    } else {
      while (erased[j]) j++;
      (*ptrs)[i] = slices->ptrs[j];
      j++;
      (*ptrs)[x] = slices->ptrs[i];
      x++;
    }
  }
  for (i = k; i < k+m; i++) {
    if (erased[i]) {
      (*ptrs)[x] = slices->ptrs[i];
      x++;
    }
  }
  return ptrs;
}

static JER_Schedule *JER_Generate_Decoding_Schedule(JER_Bitmatrix *jbm, vector <int> & erasures, int smart)
{
  int i, j, x, drive, y, index, z;
  int r,c;
  JER_Bitmatrix *decoding_matrix, *inverse, *real_decoding_matrix;
  int *row_ids;
  int *ind_to_row;
  int ddf, cdf;
  JER_Schedule *schedule;
  int *b1, *b2;
  vector <int> erased;
  int xor_result;
  int b1_row;
  int b2_row;

  int k = jbm->SUPERC;
  int w = jbm->W;
  int m = jbm->Elts.size()/w;

  if(!k || !m)return NULL;

  /* First, figure out the number of data drives that have failed, and the number of coding drives that have failed: ddf and cdf */
  ddf = 0;
  cdf = 0;
  for (i = 0; i < erasures.size(); i++) {
    if (erasures[i] < k) ddf++; else cdf++;
  }

  JER_Erasures_To_Erased(k, m, erasures, erased);
  if (erased.size() != k+m) return NULL;

  row_ids = talloc(int, k+m);
  ind_to_row = talloc(int, k+m);

  //setup id's for scheduled decoding:
  //See set_up_ptrs_for_scheduled_decoding for how these are set 
  j = k;
  x = k;
  //loop through data drives (id=0 to k-1)
  for (i = 0; i < k; i++) {
    if (erased[i] == 0) {
      //data drive i is not erased
      row_ids[i] = i;
      ind_to_row[i] = i;
    } else {
      //data drive i is erased
      while (erased[j]) j++;
      //j is now the lowest unused non-failed coding drive.
      row_ids[i] = j;
      ind_to_row[j] = i;
      j++;
      row_ids[x] = i;
      ind_to_row[i] = x;
      x++;
    }
  }
  //loop through the coding drives (id = k to k+m-1)
  for (i = k; i < k+m; i++) {
    if (erased[i]) {
      row_ids[x] = i;
      ind_to_row[i] = x;
      x++;
    }
  }

  /* Now, we're going to create one decoding matrix which is going to 
     decode everything with one call.  The hope is that the scheduler
     will do a good job.    This matrix has w*e rows, where e is the
     number of erasures (ddf+cdf) */


  real_decoding_matrix = new JER_Bitmatrix(cdf+ddf,k,w);

  /* First, if any data drives have failed, then initialize the first
     ddf*w rows of the decoding matrix from the standard decoding
     matrix inversion */
  if (ddf > 0) {

    decoding_matrix = new JER_Bitmatrix(k,k,w);

    //loop through all k super rows
    //i is the index of the super-row
    for (i = 0; i < k; i++) {
      if (row_ids[i] == i) {
        for(r=i*w;r<i*w+w;r++){
          for(c=0;c<k*w;c++){
            //setup an identity matrix in super-cell i,i
            if(r == c){
              decoding_matrix->Set(r,c,1);
            }else{
              //set all other bits in ith super row = 0
              decoding_matrix->Set(r,c,0);
            }
          }
        }
      } else {
        //the ith super-row should be the coding matrix's (row_ids[i]-k) super-row 
        //copy all w rows making up the super-row
        for(r=0;r<w;r++){
          decoding_matrix->Elts[i*w+r]=jbm->Elts[(row_ids[i]-k)*w+r];
        }
      }
    }

    decoding_matrix->Inverse();

    //set the first ddf super-rows of real_decoding matrix equal to
    //the decoding_matrix's row_ids[k+i] super-row
    for (i = 0; i < ddf; i++) {
      //w rows make up a super-row
      for(r=0;r<w;r++){
        real_decoding_matrix->Elts[i*w+r] = decoding_matrix->Elts[row_ids[k+i]*w+r];
      }
    }
    delete decoding_matrix;
  }

  /* Next, here comes the hard part.  For each coding node that needs
     to be decoded, you start by putting its rows of the distribution
     matrix into the decoding matrix.  If there were no failed data
     nodes, then you're done.  However, if there have been failed
     data nodes, then you need to modify the columns that correspond
     to the data nodes.  You do that by first zeroing them.  Then
     whereever there is a one in the distribution matrix, you XOR
     in the corresponding row from the failed data node's entry in
     the decoding matrix.  The whole process kind of makes my head
     spin, but it works.
   */
  for (x = 0; x < cdf; x++) {
    drive = row_ids[x+ddf+k]-k;

    //set real_decoding_matrix's (ddf+x)'th super-row equal to bitmatrix's drive'th super-row
    for(i=0;i<w;i++){
      //set all w rows in the super-row
      real_decoding_matrix->Elts[w*(ddf+x)+i] = jbm->Elts[w*drive+i];
    }

    //loop through all super-columns
    for (i = 0; i < k; i++) {
      if (row_ids[i] != i) {

        //set real_decoding_matrix's super-cell (ddf+x, i) = 0
        for(r=0;r<w;r++){
          for(c=0;c<w;c++){
            real_decoding_matrix->Set((ddf+x)*w+r,i*w+c,0);
          }
        }
      }
    }

    // There's the yucky part 
    //index is start of super-row drive
    //index = drive*k*w*w;
    //loop through super-columns
    for (i = 0; i < k; i++) {
      if (row_ids[i] != i) {
        //b1 is the start of super-row ind_to_row[i]-k
        //b1 = real_decoding_matrix+(ind_to_row[i]-k)*k*w*w;
        b1_row = (ind_to_row[i]-k)*w;
        for (j = 0; j < w; j++) {
          b2_row = (ddf+x)*w+j;
          for (y = 0; y < w; y++) {
            if(jbm->Isset(drive*w+j,i*w+y)){
              //loop through all the columns
              for (z = 0; z < k*w; z++) {
                xor_result = real_decoding_matrix->Isset(b2_row,z) ^ real_decoding_matrix->Isset(b1_row+y,z);
                real_decoding_matrix->Set(b2_row,z,xor_result);
              }
            }
          }
        }
      }
    }
  }

  schedule = new JER_Schedule;

  if (smart) {
    schedule = JER_BM_To_Schedule_CSHR(real_decoding_matrix);
  } else {
    schedule = JER_BM_To_Schedule_Dumb(real_decoding_matrix);
  }
  free(row_ids);
  free(ind_to_row);
  delete real_decoding_matrix;
  return schedule;
}

int JER_Schedule_Decode_Lazy(JER_Slices *slices, JER_Bitmatrix *jbm, vector <int> &erasures, int smart){

  int i, tdone;
  vector <unsigned char *> *ptrs;
  JER_Schedule * schedule;

  if(erasures.size() == 0){
	return 0;
  }

  ptrs = JER_Set_Up_Ptrs_For_Scheduled_Decoding(slices, erasures);
  if (ptrs == NULL) return -1;

  schedule = JER_Generate_Decoding_Schedule(jbm, erasures, smart);
  if (schedule == NULL) {
    delete(ptrs);
    return -1;
  }

  for (tdone = 0; tdone < slices->size; tdone += slices->PS * slices->W) {
    JER_Do_Scheduled_Operations(*ptrs, schedule, slices->PS);
    for (i = 0; i < slices->K+slices->M; i++) (*ptrs)[i] += (slices->PS * slices->W);
  }

  delete(schedule);
  delete(ptrs);

  return 0;

}

void JER_Free_Schedule_Cache(vector <JER_Schedule *> *cache){

	int e1,e2;
	int kp2; //k+2 = sqrt(cache->size())

  if(cache == NULL) return;

	//size of the cache is (k+2)^2
	kp2 = (int)sqrt(cache->size());
	if(kp2 * kp2 != cache->size()){
		return;
	}

  for (e1 = 0; e1 < kp2; e1++) {
    for (e2 = 0; e2 < e1; e2++) {
			delete (*cache)[e1*(kp2)+e2];
		}
		delete (*cache)[e1*(kp2)+e1];
  }

	delete cache;

  return;

}

//scheduling technique is 0 for dumb, 1 for CSHR
vector <JER_Schedule *> *JER_RAID_6_Schedule_Cache(JER_Bitmatrix *jbm, int technique){

  int m;
  int k;
  vector <JER_Schedule *> *scache;
  vector <int> erasures;
  int e1,e2;

  scache = new vector <JER_Schedule *>;

  m = jbm->Elts.size()/jbm->W;
  k = jbm->SUPERC;

  //only works for m = 2
  if(m != 2) return NULL;

  try{ scache->resize((k+m)*(k+m)); }
  catch(...){ return NULL; }

  /* Ok -- this is yucky, but it's how I'm doing it.  You will make an index out of erasures, which will be  e1*(k+m)+(e2).  If there is no e2, then e2 = e1.
     Isn't that clever and confusing.  Sorry.
     We're not going to worry about ordering -- in other words, the schedule for e1,e2 will be the same as e2,e1.  They will have the same pointer -- the schedule will not be duplicated. */

  erasures.resize(3,0);
  for (e1 = 0; e1 < k+m; e1++) {
      erasures.resize(2);
    erasures[0] = e1;
    for (e2 = 0; e2 < e1; e2++) {
      erasures[1] = e2;
      (*scache)[e1*(k+m)+e2] = JER_Generate_Decoding_Schedule(jbm, erasures, technique);
			(*scache)[e2*(k+m)+e1] = (*scache)[e1*(k+m)+e2];
    }
    erasures.resize(1);
    (*scache)[e1*(k+m)+e1] = JER_Generate_Decoding_Schedule(jbm, erasures, technique);
  }
  return scache;

}

int JER_Schedule_Decode_Cache(JER_Slices *slice, vector <JER_Schedule *> &cache, vector <int> &erasures){

  vector <unsigned char *>*ptrs;
  int i,k, m, w;
  int index;
  JER_Schedule * schedule;
  int size, packetsize;
  int tdone;

  k = slice->K;
  m = slice->M;
  w = slice->W;
  size = slice->size;
  packetsize = slice->PS;

  if (erasures.size() == 1){
    index = erasures[0]*(k+m) + erasures[0];
  } else if (erasures.size() == 2) {
    index = erasures[0]*(k+m) + erasures[1];
  } else {
    return -1;

  }

  schedule = cache[index];
  ptrs = JER_Set_Up_Ptrs_For_Scheduled_Decoding(slice, erasures);
  if (ptrs == NULL) return -1;

  for (tdone = 0; tdone < size; tdone += packetsize*w) {
    JER_Do_Scheduled_Operations(*ptrs, schedule,packetsize);
    for (i = 0; i < k+m; i++) (*ptrs)[i] += (packetsize*w);
  }

  delete(ptrs);

  return 0;

}

/***** End schedule functions *****/

/* Assume decoding_jm and dm_ids are already allocated just like original spec. */
int JER_Make_Decoding_Matrix(JER_Matrix *jm, vector <int> &erased, 
		JER_Matrix *decoding_jm, int *dm_ids)
{
	int i, j;
	int k, m;

	decoding_jm->W = jm->W;
	m = jm->Elts.size();
	if( m )
		k = jm->Elts[0].size();
	else
		return -1;

	/* Determine non-erased drives. Also resize decoding matrix simultaneously. */
	j = 0;
	try
	{
	decoding_jm->Elts.resize(k);
	for (i = 0; j < k; i++) {
		if (erased[i] == 0) {
			decoding_jm->Elts[j].resize( k );
			dm_ids[j] = i;
			j++;
		}
	}
	}
	catch(...) { return 1; }

	/* If the ith device is a data device, the ith row is the dm_ids[i]th row
	   of the Identity matrix. Otherwise take the dm_ids[i]-kth row of the
	   generator matrix. */
	for (i = 0; i < k; i++) {
		if (dm_ids[i] < k) {
			for (j = 0; j < k; j++) 
				decoding_jm->Elts[i][j] = 0;
			decoding_jm->Elts[i][dm_ids[i]] = 1;
		} else {
			decoding_jm->Elts[i] = jm->Elts[dm_ids[i]-k];
		}
	}

	/* Create inverse of the current decoding matrix. */
	i = decoding_jm->Inverse();
	return i;
}

// Internal Routine 
/* Assume decoding_jm and dm_ids are already allocated just like original spec. */
int JER_Make_Decoding_Bitmatrix(JER_Bitmatrix *jbm, vector <int> &erased, 
		JER_Bitmatrix *decoding_jbm, int *dm_ids)
{
	int i, j, l;
	int k, m, w;

	w = decoding_jbm->W = jbm->W;
	m = jbm->Elts.size() / w;
	k = decoding_jbm->SUPERC = jbm->SUPERC;

	/* Make sure decoding_jbm is of the right size. */
	try
	{
	decoding_jbm->Elts.resize(k*w);
	for( i = 0; i < k*w; i++ )
		decoding_jbm->Elts[i].resize( (k*w - 1) / ( 8 * sizeof( char ) ) + 1 );
	}
	catch(...){ return 1; }

	/* Determine first k surviving drives. */
	j = 0;
	for (i = 0; j < k; i++) {
		if (erased[i] == 0) {
			dm_ids[j] = i;
			j++;
		}
	}

	/* If the ith drive is a data drive, create the dm_ids[i]th super row of the
	   Identity matrix. Otherwise use the dm_ids[i]-kth super row of the generator
	   matrix. */
	for (i = 0; i < k; i++) {
		if (dm_ids[i] < k) {
			for( l = 0; l < w; l++ )
			{
				for (j = 0; j < decoding_jbm->Elts[i].size(); j++) 
					decoding_jbm->Set(i*w+l, j, 0);
				decoding_jbm->Set(i*w+l, dm_ids[i]*w+l, 1);
			}
		} else {
			for( j = 0; j < w; j++ )
				decoding_jbm->Elts[i*w+j] = jbm->Elts[(dm_ids[i]-k)*w+j];
		}
	}

	/* Invert the current decoding_jbm. */
	i = decoding_jbm->Inverse();
	return i;
}

/* Not responsible if m, k, and w in slices and jm don't match. */
int JER_Matrix_Decode(JER_Slices *slices, JER_Matrix *jm, int row_k_ones, 
		vector <int> &erasures )
{
	int i, j, edd, lastdrive;
	JER_Matrix *decoding_matrix;
	int dm_ids[slices->K];
	vector <int> erased;
	int k, m, w;

	w = slices->W;
	k = slices->K;
	m = slices->M;
	/* Only works for 8,16,32. Return with error otherwise. */
	if (w != 8 && w != 16 && w != 32) return -1;

	/* Determine erased array. */
	if( JER_Erasures_To_Erased(k, m, erasures, erased) )
		return -1;

	// Find the number of data drives failed 

	lastdrive = k;

	edd = 0;
	for (i = 0; i < k; i++) {
		if (erased[i]) {
			edd++;
			lastdrive = i;
		}
	}

	/* You only need to create the decoding matrix in the following cases:

	   1. edd > 0 and row_k_ones is false.
	   2. edd > 0 and row_k_ones is true and coding device 0 has been erased.
	   3. edd > 1

	   We're going to use lastdrive to denote when to stop decoding data.
	   At this point in the code, it is equal to the last erased data device.
	   However, if we can't use the parity row to decode it (i.e. row_k_ones=0
	   or erased[k] = 1, we're going to set it to k so that the decoding 
	   pass will decode all data.
	 */

	/* If the first coding drive is not all ones or if it was erased, lastdrive
	   equals k and we have to decode all drives with a decoding matrix. */
	if (!row_k_ones || erased[k]) lastdrive = k;

	decoding_matrix = NULL;


	/* If we have to make a decoding matrix, allocate the data first and call
	   the make_decoding_matrix method. */
	if (edd > 1 || (edd > 0 && (!row_k_ones || erased[k]))) {
		//try{ dm_ids.resize(k); }
		//catch(...){ return -1; }

		decoding_matrix = new JER_Matrix( k, k, w );
		if (decoding_matrix == NULL) {
			return -1;
		}

		if ( JER_Make_Decoding_Matrix( jm, erased, decoding_matrix, dm_ids) < 0)
		{
			free(decoding_matrix);
			return -1;
		}
	}

	/* Decode the data drives.  
	   If row_k_ones is true and coding device 0 is intact, then only decode edd-1 drives.
	   This is done by stopping at lastdrive.
	   We test whether edd > 0 so that we can exit the loop early if we're done.
	 */

	for (i = 0; edd > 0 && i < lastdrive; i++) {
		if (erased[i]) {
			JER_Matrix_Dotprod( slices, decoding_matrix, i, dm_ids, i );
			edd--;
		}
	}

	// Then if necessary, decode drive lastdrive 

	if (edd > 0) {
		for (i = 0; i < k; i++) {
			dm_ids[i] = (i < lastdrive) ? i : i+1;
		}
		JER_Matrix_Dotprod( slices, jm, 0, dm_ids, lastdrive );
	}

	// Finally, re-encode any erased coding devices 

	for( i = 0; i < k; i++ )
		dm_ids[i] = i;
	for (i = 0; i < m; i++) {
		if (erased[k+i]) {
			JER_Matrix_Dotprod( slices, jm, i, dm_ids, k+i );
		}
	}

	if (decoding_matrix != NULL) delete decoding_matrix;

	return 0;
}

/* Not responsible if m, k, and w in slices and jm don't match. */
int JER_Bitmatrix_Decode( JER_Slices *slices, JER_Bitmatrix *jbm, int row_k_ones,
		vector <int> &erasures )
{
	int i, j;
	int k, m, w;
	vector <int> erased;
	JER_Bitmatrix *decoding_matrix;
	int dm_ids[slices->K];
	int edd, *tmpids, lastdrive;

	w = slices->W;
	k = slices->K;
	m = slices->M;
	if( JER_Erasures_To_Erased(k, m, erasures, erased) )
		return -1;
	/* See jerasure_matrix_decode for the logic of this routine.  This one works 
	   just like it, but calls the bitmatrix ops instead */
	lastdrive = k;

	edd = 0;
	for (i = 0; i < k; i++) {
		if (erased[i]) {
			edd++;
			lastdrive = i;
		} 
	}

	/* If the first coding drive is not all ones or if it was erased, lastdrive
	   equals k and we have to decode all drives with a decoding matrix. */
	if (row_k_ones != 1 || erased[k]) lastdrive = k;

	decoding_matrix = NULL;

	if (edd > 1 || (edd > 0 && (row_k_ones != 1 || erased[k]))) {
	//	try{ dm_ids.resize(k); }
	//	catch(...){ return -1; }

		decoding_matrix = new JER_Bitmatrix( k, k, w );
		if (decoding_matrix == NULL) {
			return -1;
		}

		if ( JER_Make_Decoding_Bitmatrix( jbm, erased, decoding_matrix, dm_ids) < 0) {
			free(decoding_matrix);
			return -1;
		}
	}

	for (i = 0; edd > 0 && i < lastdrive; i++) {
		if (erased[i]) {
			JER_Bitmatrix_Dotprod( slices, decoding_matrix, i, dm_ids, i );
			edd--;
		}
	}

	if (edd > 0) {
		for (i = 0; i < k; i++) {
			dm_ids[i] = (i < lastdrive) ? i : i+1;
		}
		JER_Bitmatrix_Dotprod( slices, jbm, 0, dm_ids, lastdrive );
	}

	for( i = 0; i < k; i++ )
		dm_ids[i] = i;
	for (i = 0; i < m; i++) {
		if (erased[k+i]) {
			JER_Bitmatrix_Dotprod( slices, jbm, i, dm_ids, k+i );
		}
	}

	if (decoding_matrix != NULL) free(decoding_matrix);

	return 0;
}

void JER_Matrix_Encode( JER_Slices *slices, JER_Matrix *jm )
{
	int i, w;

	w = slices->W;
	if (w != 8 && w != 16 && w != 32) {
		fprintf(stderr, "ERROR: jerasure_matrix_encode() and w is not 8, 16 or 32\n");
		exit(1);
	}

	for (i = 0; i < slices->M; i++) {
		JER_Matrix_Dotprod( slices, jm, i, NULL, slices->K+i );
	}
}

void JER_Bitmatrix_Encode( JER_Slices *slices, JER_Bitmatrix *jbm )
{
	int i, w;
	int packetsize;

	packetsize = slices->PS;
	w = slices->W;

	if (packetsize%sizeof(long) != 0) {
		fprintf(stderr, "jerasure_bitmatrix_encode - packetsize(%d) %c sizeof(long) != 0\n", packetsize, '%');
		exit(1);
	}
	if (slices->size%(packetsize*w) != 0) {
		fprintf(stderr, "jerasure_bitmatrix_encode - size(%d) %c (packetsize(%d)*w(%d))) != 0\n", 
				slices->size, '%', packetsize, w);
		exit(1);
	}

	for (i = 0; i < slices->M; i++) {
		JER_Bitmatrix_Dotprod( slices, jbm, i, NULL, slices->K+i );
	}
}

/* src_ids changed to int[] instead of vector<int> so that it may be passed as 
   NULL. */
void JER_Matrix_Dotprod(JER_Slices *slices, JER_Matrix *jm, int row,
		int *src_ids, int dest_id)
{
	int init;
	unsigned char *dptr, *sptr;
	int i;
	int w, k, size;
	IVec mrow;

	mrow = jm->Elts[row];
	w = slices->W;
	k = slices->K;
	size = slices->size;

	if (w != 1 && w != 8 && w != 16 && w != 32) {
		fprintf(stderr, "ERROR: jerasure_matrix_dotprod() called and w is not 1, 8, 16 or 32\n");
		exit(1);
	}

	init = 0;

	dptr = slices->ptrs[dest_id];

	// First copy or xor any data that does not need to be multiplied by a factor 

	for (i = 0; i < k; i++) {
		if (mrow[i] == 1) {
			if( src_ids == NULL )
				sptr = slices->ptrs[i];
			else
				sptr = slices->ptrs[src_ids[i]];
			if (init == 0) {
				memcpy(dptr, sptr, size);
				jerasure_total_memcpy_bytes += size;
				init = 1;
			} else {
				galois_region_xor(sptr, dptr, dptr, size);
				jerasure_total_xor_bytes += size;
			}
		}
	}

	// Now do the data that needs to be multiplied by a factor 

	for (i = 0; i < k; i++) {
		if (mrow[i] != 0 && mrow[i] != 1) {
			if (src_ids == NULL)
				sptr = slices->ptrs[i];
			else
				sptr = slices->ptrs[src_ids[i]];
			if (src_ids == NULL) 
				sptr = slices->ptrs[i];
			else 
				sptr = slices->ptrs[src_ids[i]];
			switch (w) {
				case 8:  galois_w08_region_multiply(sptr, mrow[i], 
								 size, dptr, init); 
						 break;
				case 16: galois_w16_region_multiply(sptr, mrow[i], 
								 size, dptr, init); 
						 break;
				case 32: galois_w32_region_multiply(sptr, mrow[i], 
								 size, dptr, init); 
						 break;
			}
			jerasure_total_gf_bytes += size;
			init = 1;
		}
	}
}

/* src_ids changed to int[] instead of vector<int> so that it may be passed as 
   NULL. */
void JER_Bitmatrix_Dotprod( JER_Slices *slices, JER_Bitmatrix *jbm, 
		int srow, int *src_ids, int dest_id )
{
	unsigned char *bpptr;
	int size, w, packetsize, k;
	char first;
	int i, j, sindex;

	size = slices->size;
	packetsize = slices->PS;
	w = slices->W;
	k = slices->K;

	if (size%(w*packetsize) != 0) {
		fprintf(stderr, "jerasure_bitmatrix_dotprod - size%c(w*packetsize)) must = 0\n", '%');
		exit(1);
	}

	bpptr =  slices->ptrs[dest_id];

	// For each set of packets of Dt
	for (sindex = 0; sindex < size; sindex += (packetsize*w)) {
		// For each packet in the current set of Dt
		for( i = 0; i < w; i++ )
		{
			first = 1;
			for( j = 0; j < k * w; j++ )
			{
				if( jbm->Isset(srow*w+i, j) )
				{
					if( first )
					{
						if( src_ids == NULL )
							memcpy( bpptr + i*packetsize + sindex,
									slices->ptrs[j/w]+
									(j%w)*packetsize+sindex, packetsize );
						else
							memcpy( bpptr + i*packetsize + sindex,
									slices->ptrs[src_ids[j/w]]+
									(j%w)*packetsize+sindex, packetsize );
						jerasure_total_memcpy_bytes += packetsize;
						first = 0;
					}
					else
					{
						if( src_ids == NULL )
							galois_region_xor( bpptr + i*packetsize + sindex,
									slices->ptrs[j/w]+
									(j%w)*packetsize+sindex, 
									bpptr + i*packetsize + sindex,
									packetsize );
						else
							galois_region_xor( bpptr + i*packetsize + sindex,
									slices->ptrs[src_ids[j/w]]+
									(j%w)*packetsize+sindex, 
									bpptr + i*packetsize + sindex,
									packetsize );
						jerasure_total_xor_bytes += packetsize;
					}
				}
			}
		}
	}
}

void JER_Do_Parity( JER_Slices *slices, unsigned char *parity_ptr )
{
	int i;
	int size, k;
	
	size = slices->size;
	k = slices->K;

	memcpy(parity_ptr, slices->ptrs[0], size);
	jerasure_total_memcpy_bytes += size;

	for (i = 1; i < k; i++) {
		galois_region_xor(slices->ptrs[i], parity_ptr, parity_ptr, size);
		jerasure_total_xor_bytes += size;
	}
}

/* Converts a list-style version of the erasures into an array of k+m elements
   where the element = 1 if the index has been erased, and zero otherwise */
int JER_Erasures_To_Erased(int k, int m, vector <int> &erasures, vector <int> &erased)
{
	int td;
	int t_non_erased;
	int i;

	try{ erased.resize( k + m ); }
	catch(...){ return -1; }
	t_non_erased = k + m;

	for (i = 0; i < k + m; i++) erased[i] = 0;

	for (i = 0; i < erasures.size(); i++) {
		if (erased[erasures[i]] == 0) {
			erased[erasures[i]] = 1;
			t_non_erased--;
			if (t_non_erased < k) {
				return -1;
			}
		}
	}
	return 0;
}

JER_Matrix* JER_Matrix_Multiply( JER_Matrix* m1, JER_Matrix* m2 )
{
	int i, j, k, l;
	JER_Matrix* product;

	int r1, r2, c2, w;

	r1 = m1->Elts.size();
	if( r1 <= 0 || m1->Elts[0].size() != m2->Elts.size() || 
		m2->Elts.size() <= 0 )
		return NULL;
	c2 = m2->Elts[0].size();
	r2 = m2->Elts.size();
	w = m1->W;
	product = new JER_Matrix( r1, c2, w );
	if( product == NULL )
		return NULL;

	for (i = 0; i < r1; i++) {
		for (j = 0; j < c2; j++) {
			for (k = 0; k < r2; k++) {
				product->Elts[i][j] ^= galois_single_multiply(m1->Elts[i][k], 
															  m2->Elts[k][j], w);
			}
		}
	}
	return product;
}

void JER_Get_Stats( vector<double> &fill_in)
{
	try	{ fill_in.resize(3); }
	catch(...){	return; }

	fill_in[0] = jerasure_total_xor_bytes;
	fill_in[1] = jerasure_total_gf_bytes;
	fill_in[2] = jerasure_total_memcpy_bytes;
	jerasure_total_xor_bytes = 0;
	jerasure_total_gf_bytes = 0;
	jerasure_total_memcpy_bytes = 0;
}
